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ABSTRACT 

Motivated by suggestions that binaries with almost equal- mass components ("twins") play an im- 
portant role in the formation of double neutron stars and may be rather abundant among binaries, 
we study the stability of synchronized close and contact binaries with identical components in circular 
orbits. In particular, we investigate the dependency of the innermost stable circular orbit on the core 
mass, and we study the coalescence of the binary that occurs at smaller separations. For twin binaries 
composed of convective main-sequence stars, subgiants, or giants with low mass cores (Mc < 0.15M, 
where M is the mass of a component), a secular instability is reached during the contact phase, accom- 
panied by a dynamical mass transfer instability at the same or at a slightly smaller orbital separation. 
Binaries that come inside this instability limit transfer mass gradually from one component to the 
other and then coalesce quickly as mass is lost through the outer Lagrangian points. For twin giant 
binaries with moderate to massive cores {Mc > 0.15Af), we find that stable contact configurations 
exist at all separations down to the Roche limit, when mass shedding through the outer Lagrangian 
points triggers a coalescence of the envelopes and leaves the cores orbiting in a central tight binary. 
In addition to the formation of binary neutron stars, we also discuss the implications of our results 
for the production of planetary nebulae with double degenerate central binaries. 
Subject headings: binaries: close — binaries: general — hydrodynamics — instabilities — methods: 
numerical — stars: general 



I. INTRODUCTION AND MOTIVATION 
1.1. Formation of Binary Neutron Stars 

The evolutionary history and formation of close bina- 
ries with two neutron stars similar to the Hulse- Taylor 
pulsar B1913-I-16 (Hulse & Taylor 1975) and the dou- 
ble pulsar J0737-3039 (Burgay et al. 2003) is a topic 
of intense current interest. Most recent studies of the 
known double neutron stars focus on the stages going 
back to the time of the second supernova explosion and 
the formation of the youngest of the two neutron stars 
(Dcwi & van den Hcuvcl 2004; Willems & Kalogera 2004; 
Willems et al. 2004; Piran & Shaviv 2005; Stairs et al. 
2006; Wang, Lai, & Han 2006; Wong et al. 2010, and 
references therein). Although these studies provide very 
interesting constraints on the properties of the stellar 
progenitor of the second neutron star, they do not probe 
the earlier evolutionary history. That part remains un- 
certain and more difficult to constrain empirically based 
on the measured properties of observed systems. 

Since the discovery of the Hulse- Taylor binary, the ori- 
gin of double neutron stars has been naturally connected 
to the evolution of massive binaries, with stellar com- 
ponents that are massive enough to form two neutron 
stars at the end of their lifetime. Over the years, a qual- 
itative consensus of understanding for the formation of 
double neutron stars developed (see, e.g., Bhattacharya 
& van den Heuvel 1991): massive binaries experience a 
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phase of stable mass transfer when the primary over- 
flows its Roche lobe revealing its helium core; this core 
ends its lifetime in a supernova forming the first neu- 
tron star in the system; the binary becomes a high-mass 
X-ray binary until the massive secondary fills its Roche 
lobe and the binary enters a dynamically unstable phase 
of mass transfer leading to inspiral in a common enve- 
lope that engulfs the neutron star; during this phase the 
neutron star is thought to be spun up through recycling 
and, if the binary avoids a merger in the inspiral, the 
helium core of the secondary is revealed after the enve- 
lope ejection; this core explodes in a supernova, forming 
the second neutron star in the system, and the double 
neutron star further evolves through orbital contraction 
and gravitational-wave emission. Variations of this evo- 
lutionary sequence have been shown to be realized by 
theoretical binary population studies (e.g., Belczynski et 
al. 2002). 

Brown (1995) argued that the inspiral of the neutron 
star during the common envelope phase in the standard 
model is problematic: the neutron star is expected to ex- 
perience hypercritical accretion (Chevalier 1993) at rates 
many orders of magnitude above the photon Eddington 
limit through a neutrino-cooled accretion flow (see also 
Fryer et al. 1996). Such a rapid accretion phase along 
with the adoption of a low maximum neutron star mass 
(~1.5M0 derived for a soft equation of state for neutron 
star matter with kaon condensation) led Brown to con- 
clude that all neutron stars in common envelope phases 
will accrete enough matter to collapse into a black hole. 
Consequently, he argued that the "standard" evolution- 
ary sequence described above aborts the formation of 
double neutron stars and instead leads to the formation 
of binaries with low-mass black holes and neutron star 
companions. We note however that, even with the same 
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treatment of hypercritical accretion, a maximum neu- 
tron star mass of ^ 2Mq (corresponding to a more reg- 
ular stiff equation of state) does prevent a good fraction 
of neutron stars from being transformed into low-mass 
black holes (Belczynski et al. 2002). 

Brown (1995) also noted how the masses measured in 
known double neutron stars are very close to being equal 
(Nice et al. 1996; Thorsett & Chakrabarty 1999; Stairs 
et al. 2002; Weisberg & Taylor 2005; Jacoby et al. 2006; 
Kramer et al. 2006, and references therein). Motivated 
by these two points, he proposed a different formation 
channel for double neutron stars. Brown suggested that 
double neutron stars form from massive binaries with 
component masses that are within ^4% of one another. 
Consequently the red giant phases of the two compo- 
nents coincide in time and when mass transfer ensues 
from the primary, both components have deep convective 
envelopes and well developed helium cores, so that a dou- 
ble core phase develops, where the two helium cores orbit 
within the combined envelopes of the two massive stars. 
Provided that there is enough orbital energy, the com- 
mon envelope is ejected before the two cores merge, and 
a tight binary with two helium cores is formed. These 
two cores differ very little in mass and reach core col- 
lapse one very soon after the other (~ 10^ yr based on 
helium-star models, Habets 1986; Pols 1994), forming a 
close double neutron star. 

The advantages of this hypothesized evolutionary 
channel are (i) a neutron star never experiences com- 
mon envelope spiral-in and hypercritical accretion, and 
(ii) the two stellar components have so similar masses 
that they naturally form neutron stars of almost equal 
mass as observed (Bethe & Brown 1998; Bethe, Brown, 
& Lee 2007). On the other hand, this channel (i) re- 
quires that mass transfer between the red giant progeni- 
tors will indeed lead to the inspiral of the two cores in a 
common envelope and (ii) requires fine-tuning the condi- 
tions for recycling the first neutron star, as this spin up 
must occur during the very short interval between the 
two supernovae through the stellar wind or possibly dur- 
ing the brief Roche-lobe overflow from the lower-mass 
helium star (Dewi et al. 2006). A potential additional 
disadvantage is that this channel is very restrictive in 
that it requires progenitors that are at most only ^ 4% 
apart in mass: however, in their study of protobinary 
stars, Krumholz & Thompson (2007) find that for a wide 
range of initial conditions Roche lobe overflow tends to 
equalize the masses of the binary components. 

The double neutron star formation channel suggested 
by Brown (1995) has attracted renewed attention be- 
cause of the reported abundance of "twins," massive bi- 
naries with mass ratios very close to unity (within 5%) 
by Pinsonneault & Stanek (2006). Specifically they ana- 
lyze data for 21 detached eclipsing binaries in the Small 
Magellanic Cloud and find that the data are consistent 
with a flat mass hmction containing 55% of the systems 
and "twins" with mass ratios greater than 0.95 contain- 
ing the remaining 45% of the population. However, it is 
important to note that there are severe selection effects 
against the discovery of binaries with small mass ratios 
(typically < 0.5, Hogeveen 1992a, b,c); therefore the con- 
tribution of twins may not be as significant as implied by 
the most recent observations. Quantitative modeling of 
the associated selection effects is required to derive more 



reliable statistical conclusions. 

Apart from uncertainties with the initial properties of 
the binary population, a number of physical processes 
related to these formation channels make it hard to as- 
sess their relative contribution to double neutron star 
formation. The physics of common envelopes and hyper- 
critical accretion, as well as of neutron star equations of 
state and their maximum mass, is not well understood 
(but see Lee, Park, & Brown 2007). Also, the develop- 
ment of a dynamical instability and subsequent common 
envelope phase with the inspiral of the two cores is as- 
sumed, but has not been investigated before in any de- 
tail. In this study, we attempt to imderstand better one 
of the aspects related to the formation chaniic^l suggested 
by Brown (1995): the fate of mass transfer between two 
stars of almost equal mass. 

1.2. Binaries and Planetary Nebulae 

Although our primary motivation is to investigate 
a formation channel for binary neutron stars, our re- 
sults are relevant to other scenarios as well. For exam- 
ple, through common envelope evolution, binaries could 
transform quite naturally into planetary nebulae (PNe) 
with a central close binary. The influence of binaries on 
the production and morphology of PNe has received in- 
creased attention in recent years (e.g., Miszalski et al. 
2009; Jones et al. 2010), and observations to date have 
identified approximately 40 close binaries in the centers 
of PNe: see De Marco (2009) for a summary of both the 
relevant theory and observations. As lifetimes of PNe 
are less than only 10"'' years, an embedded central binary 
has undergone no significant evolution since the common 
envelope phase that presumably formed it. 

About a quarter of the known close binary systems 
within PNe are thought to be double degenerates, that 
is, binaries in which the components are pre-white dwarfs 
(also known as extreme horizontal branch stars, hot sub- 
dwarfs, or subdwarf B and O stars) or white dwarfs. 
Such double degenerate binaries, with or without plan- 
etary nebulae, are particularly interesting. As possible 
progenitors for Type la supernovae, they are the desired 
targets of the ESO SN la Progenitor Survey (see Geier 
et al. 2011, and references therein). 

Interestingly, three of the five well studied double de- 
generate binaries within PNe (see Hillwig 2011, and ref- 
erences therein) may contain nearly equal mass compo- 
nents and therefore are possibly the progeny of twin bi- 
naries. In particular, models of a double hot subdwarf 
binary with g « 1 are consistent with photometric and 
spectroscopic observations of the central stars in NGC 
6026 (Hillwig et al. 2010) and in Abell 41 (Bruch et al. 
2001; Shimanskii et al. 2008). Similarly, recent observa- 
tions of the central star in Hen 2-428 reveal it also to be a 
double degenerate binary, with the effective temperature 
of the components being within a few thousand Kelvin 
(Santander-Garcia et al. 2011), suggesting a mass ratio 
near g' = 1. As we investigate in this paper, the coales- 
cence of a twin giant binary could lead to the formation 
of double degenerate core surrounded by a circumbinary 
envelope, a type of proto-PN. 

1.3. Theoretical Work on Close Binaries 

Most of the classical theoretical work on close bina- 
ries was done in the limit of a self-gravitating incom- 
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pressiblo fluid (see Chandrasekhar 1969, and references 
therein). An essential result found in the incompress- 
ible case is that the hydrostatic equilibrium solutions for 
sufRciently close binaries can become globally unstable 
(Chandrasekhar 1975; Tassoul 1975). The classical an- 
alytic studies for binaries were extended to polytropes 
in the work of Lai et al. (1993a, b, 1994a, b,c). In their 
approach, the stars are modeled as self-gravitating com- 
pressible ellipsoids, and an energy variational principle 
is used to construct approximate equilibrium configura- 
tions and study their stability. These treatments, along 
with complementary numerical hydrodynamic calcula- 
tions (Rasio & Shapiro 1992, 1994, 1995), demonstrate 
that dynamical instabilities persist in the compressible 
regime and can cause a binary to coalesce to form a 
rapidly rotating spheroidal object. Such a dynamical in- 
stability can trigger a merger in just a few orbital periods 
(Rasio & Shapiro 1992) or an episode of mass transfer 
that lasts many orbits (Motl et al. 2002; D'Souza ct al. 
2006; Frank 2008; Dan et al. 2009; Loren-Aguilar et al. 
2009). 

The evolution of a close binary system can also be af- 
fected by another type of global instability. It has been 
referred to by various names, such as the secular instabil- 
ity (Lai et al. 1993a, b, 1994a, b,c), tidal instability (Coun- 
selman 1973; Hut 1980), gravogyro instability (Hachisu 
& Eriguchi 1984), and Darwin instability (Levine et al. 
1993). Its physical origin is easy to understand (Lai et al. 
1993a, 1994b; Rasio 1994; Webbink 2006). There exists 
a minimum value of the total angular momentum J for 
a synchronized close binary. This is simply becaiise the 
spin angular momentum, which increases as the separa- 
tion r decreases for a synchronized system, can become 
comparable to the orbital angular momentum for suffi- 
ciently small r. A system that reaches the minimum of 
J cannot evolve further by angular momentum loss and 
remain synchronized. Instead, the combined action of 
tidal forces and viscous dissipation will drive the system 
out of synchronization and cause rapid orbital decay as 
angular momentum is continually transferred from the 
orbit to the spins. The orbital decay then proceeds on 
a timescale comparable to the synchronization time of a 
stable binary. 

In this paper, we pay particular attention to the onset 
of orbital instabilities, including those characterized by 
mass transfer, as well as the subsequent inspiral of the 
cores. We extend previous hydrodynamic studies of close 
binary systems to cases that involve identical giants, that 
is, to twin stars with dense stellar cores and extended 
envelopes. So that our results can be scaled to stars of 
any size or mass, we approximate each star as a con- 
densed polytrope, namely, as a point mass surrounded 
by a uniform specific entropy fluid of adiabatic index 
r — 5/3. Such a model is appropriate for a fully convec- 
tive monatomic ideal gas surrounding a compact core. 
We consider fractional core masses rric that cover the 
entire range of theoretical possibilities. Zero core mass 
models correspond to normal (or "complete") n = 1.5 
polytropes, which are most appropriate for low mass 
main-sequence stars and non-relativistic white dwarfs. 
At the other extreme, when nic = 1, there is no mass 
in the gaseous envelope, and the binary is simply that 
of two point masses. By varying the core mass between 
these extremes, we are able to study in a systematic way 



the full parameter space of twin binaries and determine 
under what conditions mass transfer develops or an in- 
nermost stable circular orbit exists. 

For typical compositions, the subgiant phase begins 
when the core mass grows to ^ 10% of the total mass, 
the so-called Schonberg-Chandrasekhar limit (Schonberg 
& Chandrasekhar 1942). By the time the star reaches the 
base of the red giant branch, the core mass has increased 
by a few more percent of the total mass. Roughly speak- 
ing then, our models with nic < 0.1, 0.1 ^ rric < 0.13, 
and nic ^ 0.13 correspond to main sequence, subgiant, 
and giant stars, respectively. The assumption that the 
stellar envelope has constant specific entropy makes our 
models most relevant to stars that are fully convective 
(as with low mass main sequence stars) or that have deep 
convective envelopes (as in many red giants). 

The use of condensed polytropes as a model of red 
giants has a rich history, including seminal work by 
Chandrasekhar (1939), Osterbrock (1953), and Harm & 
Schwarzschild (1955). Mass transfer in close binary sys- 
tems has been modeled with the help of condensed poly- 
tropes as well: the response of the donor due to mass loss 
is considered, with its resulting contraction, or expan- 
sion, compared against that of its Roche lobe (Paczyihski 
& Sienkiewicz 1972; Hjellming & Webbink 1987).' If at 
the onset of mass transfer the star contracts less rapidly 
than its Roche lobe (or expands more rapidly than it), 
then the mass exchange is dynamically unstable, and 
the binary will evolve rapidly toward a new, often qual- 
itatively different, equilibrium. Hjellming & Webbink 
(1987) find that equal mass condensed polytrope con- 
tact binaries with fractional core masses rric ^ 0.46 ex- 
perience stable mass transfer. More recently, Krumholz 
& Thompson (2007) have used condensed polytropes to 
study the formation of twin star systems. 

Although essential for a qualitative understanding of 
mass transfer, such treatments of close binaries do, how- 
ever, make several simplifying approximations: most im- 
portantly, (i) the dynamics of the orbit and size of the 
Roche lobe are treated in the point mass approximation, 
(ii) the response of the binary components to mass loss 
or gain is modeled as if each star were spherical and in 
isolation, and (iii) mass that overflows a Roche lobe is 
considered to leave that star. These approximations are 
quite reasonable for semidetached binaries, but their va- 
lidity can be questioned for contact binaries. In such 
cases, a common envelope persists in equilibrium out- 
side of the Roche lobes (the inner Lagrangian surface) 
so that the pressure and density on the Roche lobes are 
non-zero: mass that overflows a Roche; lobe is not neces- 
sarily transferred to the other star but rather can persist 
in equilibrium inside the outer Lagrangian surface. A 
primary goal of this paper is therefore to relax the ap- 
proximations of previous works by using accurate hydro- 
dynamical calculations to study contact binary systems. 

Our paper is organized as follows. In §2 we review 
our numerical method and general conventions. In §3 we 
present our results for the equilibrium and stability prop- 
erties of twin binary systems. The dynamical evolution 
to complete coalescence is followed for several unstable 
systems. Implications of our results are discussed in §4. 

2. NUMERICAL METHODS AND ASSUMPTIONS 
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2.1. The SPH code 

To generate our models, we use a modified version of 
the SPH code originally developed by Rasio (1991) that 
has been updated to include the variational equations 
of motion derived in Gaburov ct al. (2010). SPH is a 
Lagrangian particle method, meaning that the fluid is 
represented by a finite number of fluid elements or "par- 
ticles." Associated with each particle i are, for example, 
its position rj, velocity v^, and mass mj. Each particle 
also carries a purely numerical smoothing length hi that 
determines the local spatial resolution and is used in the 
calculation of fluid properties such as acceleration and 
density. Details of our SPH code, such as the particular 
form of the artificial viscosity Tiij and smoothing kernel 
Wij implemented, are described in Gaburov et al. (2010). 
See Rasio & Lombardi (1999) and Rosswog (2009) for re- 
views of SPH. 

Because the gas in om initial stellar models is of con- 
stant specific entropy, wc find it convenient to integrate 
the so-called entropic variable Ai of each particle i. The 
entropic variable is simply the proportionality constant 
in the polytropic equation of state p = Ap^ ^ where p 
is pressure and p is density. The entropic variable is 
so named because of its close connection to entropy: 
both quantities are conserved in reversible processes and 
strictly increase otherwise. We therefore use dAi/dt = 
in the relaxations of our single star models and in the 
calculations of our binary equilibrium sequences. For 
our dynamical calculations of merger scenarios (see §2.4 
and §3.2), we evolve A^ according to the discretized SPH 
version of the first law of thermodynamics: 

^ = ^"T^ mjUij (v, - ) • ViWij {hi) . (1) 

Hi j 

To calculate the gravitational accelerations and po- 
tentials, we use direct summation on NVIDIA graphics 
cards, softening with the usual SPH kernel as in Hern- 
quist & Katz (1989). The use of such a softening with 
finite extent (as opposed, for example, to Plummer soft- 
ening) increases the accuracy and stability of our SPH 
models, consistent with the studies of Athanassoula et 
al. (2000) and Dehnen (2001). The gravity of core points 
in our models is similarly softened, applying a constant 
smoothing length comparable to the minimum smooth- 
ing length in the system. 

2.2. Single Star Models 

In this section we present our procedure for modeling 
the stars that are used in the binary simulations of §3. 
Hydrodynamically, a subgiant or giant can be treated 

as a two-component system: a high-density, degenerate 
core, surrounded by an extended envelope. The very 
large density contrast between the core and the enve- 
lope, along with the small core radius, justifies the use 
of a single point mass to represent the core. Because 
the giants we wish to model mostly have deep convec- 
tive envelopes with an equation of state dominated by 
monatomic ideal gas pressure, we treat their g as as a 
constant specific entropy fiuid with an adiabatic index 

r = 5/3. 

Specifically, our stellar models are the so-called con- 
densed polytropes, namely, constant specific entropy 



fluid surrounding a point mass core (Chandrasekhar 
1939; Harm & Schwarzschild 1955), which we parame- 
terize by the core mass TOc- Each of the condensed poly- 
tropes in our family of models has total mass M — 1 
and radius R = 1. The unit system is completed by 
choosing Newton's gravitational constant G = 1. Con- 
densed polytropc models have not only the advantage of 
reproducibility but also of scalability to any stellar mass 
and radius: although our focus here is on stars massive 
enough ultimately to yield a neutron star, the same cal- 
culations are also valid for low mass systems. In the 
limiting scenarios of nic = and rric = 1, we recover the 
well-studied cases of a n = 3/2 polytrope and a point 
mass, respectively. 

Figure 1 shows the pressure and density profiles as a 
function of radius for condensed polytropes with core 
masses 0, 0.125, and 0.5. For comparison, we also dis- 
play the profiles of red giant stars computed using the 
TWIN stellar evolution code (Eggleton 1971; Glebbeek, 
Pols, & Hurley 2008) from the MUSE software environ- 
ment^ (Portegies Zwart et al. 2009): we evolve 10 and 
25 Mq stars with initial helium abundance Y = 0.28 
and metallicity Z = 0.02 until they obtain core masses 
of approximately nic = 0.16 and 0.27, respectively. This 
corresponds to our initially 10 Mq star being at the very 
tip of the red giant branch and the initially 25 Mq star 
being on the upper part of the branch. The profiles of the 
simple condensed polytrope models follow the same gen- 
eral trend as those of the red giants, indicating that our 
simple models can indeed help provide an understanding 
of real red giant binaries. 

Table 1 presents the models used in this paper. Col- 
umn (1) gives the core mass rric, while column (2) lists the 
corresponding Eq value: this is the parameter E in the 
original notation of Osterbrock (1953), as well as in the 
notation of Harm & Schwarzschild (1955) and Hjcllming 
& Webbink (1987). (We prefer to reserve the variable E 
for energy.) The value of Eq controls the shape of the 
density profile: to lowest order near the surface (r « i?), 
p(r) « (2/5)3/2£;o(l - r/i?)3/2M/(47ri?3). In practice, 
we adjust Eq to achieve the desired core mass rric. 

TABLE 1 
Parent Star Characteristics 



nic 


Eo 


(1) 


(2) 





45.4808 


0.05 


39.9250 


0.1 


35.2403 


0.125 


33.1661 


0.15 


31.2407 


0.175 


29.4461 


0.2 


27.7673 


0.25 


24.7072 


0.3 


21.9782 


0.4 


17.2861 


0.5 


13.3582 


0.6 


9.9904 


0.7 


7.0496 


0.8 


4.4443 


0.9 


2.1092 


0.99 


0.1984 



We begin by making an SPH model of a single star 
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Fig. 1. — The solid curves show the quasi-analytic radial profiles 
of the pressure p and density p for nic =0,0.125, and . 5 condensed 
polytropes; for comparison, the dashed and dotted curves represent 
red giants, as modeled by the TWIN stellar evolution code, with 
initial masses respectively of 10 and 25 A^q. For the condensed 
polytropes, curves associated with larger nic are higher on the left 
edge of the figure and lower on the right edge. Units are such that 
G = M = R=l. 

in isolation. Unless stated otherwise, we use N = 
19938 SPH particles initially placed on a hexagonal close 
packed lattice with a lattice spacing constant ai — 
0.0542, with particles extending out to a radius that is 
between one and two smoothing lengths less than the full 
stellar radius. Wc model the stellar core as a point mass 
that interacts gravitationally, but not hydrodynamically, 
with the rest of the system, as suggested by Rasio & 
Shapiro (1991) and others. The gravitational influence 
of these core points are softened according to the SPH 
kernel with h = 0.0498. Particle masses are first appor- 
tioned according to the desired density profile and then 
slightly rescaled before relaxation begins to ensure that 
the correct total mass M = 1 is precisely achieved. The 
entropic variable Ai of each SPH particle is set to the 
desired polytropic constant K, which is determined from 
the mass-radius relation for n = 3/2 condensed poly- 
tropes (Osterbrock 1953): 

R = {An)-y^G-^KE^^^M-'/\ (2) 

After the initial parameters of the particles have been 
assigned, we relax the SPH fluid into hydrostatic equilib- 
rium. This relaxation is effected by including an artificial 
viscosity contribution in the acceleration equation (with 
a — 1 and /3 = 2 in equation (A19) of Gaburov et al. 
2010), while still keeping Ai constant for all particles. In 
this way, the entropy of the system is preserved while the 
system approaches equilibrium. The total energy typi- 
cally decreases by less than a percent in the process, in- 
dicating that our initial assignment of particle properties 
was indeed very close to an equilibrium state. 

This approach allows us to model the desired profiles 



Fig. 2. — Properties of the model with core mass rric = 0.25 as a 
function of radius, after relaxation for 500 time units. The frames 
in the left column compare calculated pressure p and density p pro- 
files of the star (dashed curves) against particle data from our SPH 
model (dots). The right column provides additional SPH particle 
data: individual SPH particle mass m, smoothing length h, num- 
ber of neighbors A'^jv, and radial component of the hydrodynamic 
acceleration Ohydro (upper data) and gravitational acceleration g 
(lower data). 

very accurately. An example is presented in Figure 2, 

where we plot desired profiles and SPH particle data for 
our relaxed nic = 0.25 star. Although the core and sur- 
face of the star of course cannot be resolved on the length 
scale of a smoothing length (typically 0.05 to 0.08 length 
units), the thermodynamic profiles of the SPH model 
nicely reproduce the quasi-analytic curves. Indeed, the 
SPH data in the left column of Figure 2 are difficult to 
distinguish from the desired pressure and density pro- 
files throughout most of the star. We also note that the 
hydrodynamic and gravitational accelerations are very 
nearly eqtial in magnitude and opposite in direction, as 
necessary for hydrostatic equilibrium. 

2.3. Binary Equilibrium Configurations 
2.3.1. SPH Calculations 

The ability of our code to model close and contact bi- 
nary systems for thousands of orbits or longer is pre- 
sented in Gaburov et al. (2010). Here we present our 
methods for modeling equilibrium sequences of twin bi- 
naries, that is, binaries that consist of two identical stars 
in synchronized orbit. First, we place identical relaxed 
stellar models along the x axis with their centers of mass 
separated by r. While the relaxation of the binary takes 
place, the entropic variable particle values are held con- 
stant. The center of mass of the entire system is fixed at 
the origin. In addition, the positions of the particles are 
continuously adjusted (by a simple uniform translation 
along the binary axis), so that the separation between 
the centers of mass equals the desired separation r. 

The orbit is chosen to occur in the xy plane. The an- 
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gular velocity florb defining the corotating frame is up- 
dated at every tiniestep so that the centrifugal and iner- 
tial accelerations acting on the fluid cancel. In particu- 
lar, we wish to find configurations in which r2^j.^,(a;jX + 
Uiy) = —{vx,i^ + Vy^if), where the velocity derivatives 
on the right hand side are components of acceleration 
in the inertial frame. By taking the dot product with 
mi{xiX + Uiy) and then summing over all particles, we 
obtain 

~Y.i'mi{xiVx,i+ ViVy,i) 



orb 



We also include a drag force that opposes the veloc- 
ity and provides a contribution to the acceleration of 
— Vi/troiax- We use trciax — 3, approximately the funda- 
mental period of oscillation for our parent models. We 
do not include any artificial viscosity contribution when 
finding binary equilibrium configurations. 

In order to find c^qTiilibrium configurations for a pre- 
cisely equal mass binary, even for configurations unsta- 
ble to mass transfer, we enforce a symmetry in particle 
properties: for each particle i in star 1, there is a part- 
ner particle j in star 2 at Xj = —Xi and yj = —yi with 
velocity components Vxj = —Vx,i, Vyj = —Vyj and with 
acceleration components Vxj — —Vx,i, Vyj — —Vy,i- All 
other properties are identical for any such pair of parti- 
cles. 

The separation r between the centers of mass can be 
allowed to drift slowly so that an equilibrium sequence 
is constructed: a so-called "scanning run." In practice, 
we start runs that will scan over separations by hold- 
ing the centers of mass fixed at an initial separation 
r(0) for 40 time units, allowing the system to approach 
a tidally bulged equilibrium configuration. At an addi- 
tional amount of time t, the separation is set according 

to rit) = r(0) [r(t,can)/r(0)]*/*==°" . This form for r{t) al- 
lows the change in r to occur at a decreasing rate as the 
stars approach and interact more strongly, although the 
exact form is not critical to our results. We typically use 
tscan = 300, r(0) = 3.3, and r{tscan) = 2.1. 

2.3.2. Data Reduction Methods 

Once an SPH binary calculation has completed, we 
analyze the system at various separations along the se- 
quence. To this end, a useful quantity to consider is the 
effective potential, calculated as 

$e(ar, y, z) = $(ar, y, z) - ^nl^^{x'^ + y^), (4) 

where $ is the gravitational potential, the coordinate 
y measures perpendicular to the binary axis in the or- 
bital plane, and z measures parallel to the rotation axis. 
Along the binary axis {y = z = Q), the effective potential 

has a local maximum <bf' at x = (the inner Lagrangian 

point) and global maxima $1°^ at = Xq (the outer La- 
grangian points). There are two minima at \x\ = Xc, cor- 
responding to the cores of the two components.^ In equi- 
librium, the fluid will fill up the effective potential well 
to some maximum, constant level Borrowing the 

® Note that in general r 7^ 2xc, because we define the binary 
separation r as the distance between the centers of mass of the 
two components. 



terminology from models of W UMa binaries (Rucinski 
1992, and references therein), we follow Rasio & Shapiro 
(1995) and define the degree of contact 77 as 



(5) 



Clearly, we have rj < for detached configurations: that 
is, none of the fluid has a large enough effective potential 
energy to exceed the effective potential energy barrier at 
the inner Lagrangian point. For < < 1, the effec- 
tive potential of the fiuid near a; = does exceed the 
barrier, and the system is classified as a contact binary. 
For ry > 1, the envelopes overfiow beyond the outer La- 
grangian surface, and no dynamical equilibrium configu- 
ration can be achieved; that is, the system has exceeded 
the Roche limit. Calculations in which we slowly scan 
to smaller separations can therefore determine position 
of first contact (77 = 0), the secular stability limit (at 
the minimum energy and angular momentum along the 
sequence), and the Roche limit (77 — 1). 

It is important to realize that the equilibrium sequence 
of a twin binary passes smoothly from detached to con- 
tact configurations as the separation r decreases. This 
is in contrast to all binary equilibrium sequences with 
mass ratio q ^ 1, which terminate at a Roche limit cor- 
responding to the onset of mass transfer through the in- 
ner Lagrangian point (that is, once one of the binary 
components overflows its Roche lobe). For twin binaries, 
however, the Roche limit, which we still define as the last 
equilibrium configuration along a sequence with decreas- 
ing r, corresponds to the onset of mass shedding through 
the outer Lagrangian points: as an example, note that 
several particles have been shed to the far left and far 
right of the r = 2.22 frame in Fig. 3. 

We estimate ^i^^ from our SPH models by finding the 
maximum effective potential of the points along the x- 
axis that are within one smoothing length of the cen- 
ter of an SPH particle. Thus, even if the centers of all 
SPH particles are within the outer Lagrangian surface, 
we may still consider the system as having reached the 
Roche limit when some smoothing kernels extend sub- 
stantially beyond the outer Lagrangian surface. Such an 
estimate accounts for the fact that an SPH particle is 
not a point mass but instead represents a parcel of fluid 
with a density profile described by the smoothing kernel. 

Our means of estimating $1"^ allow our critical separa- 
tion results to converge quickly to a steady value as the 
resolution is increased up to the resolution presented in 
this paper (see §3.1). 

2.4. Dynamical Calculations 

We generate initial conditions for our dynamical runs 
by taking a configuration at the desired separation r from 
a scanning run and then relaxing for an additional 200 
time units. If a particle escapes past an outer Lagrangian 
point during this time interval, then we end the relax- 
ation stage and begin following the dynamics immedi- 
ately. During dynamical calculations, we include no drag 
force and move the particles according to their velocities 
in the usual way (for details, see Gaburov et al. 2010). 
Particles are again treated independently so that mass 
transfer events can be followed; that is, unlike the scans 



Twin Binaries 



7 



described in §2.3, no symmetry constraints are applied 
to particle motion. Artificial viscosity is implemented 
in both the acceleration and the entropy equations. To 
minimize the spurious effects of artificial viscosity (Lom- 
bard! et al. 1999), our dynamical calculations are done 
in a rotating frame, with the the angular velocity fiorb 
calculated once at the beginning of the dynamical evolu- 
tion and thereafter held constant when applying Coriolis 
and centrifugal forces. 

3. RESULTS 

Using the methods described above, we construct twin 
binary sequences for 16 different core masses rric listed in 
Table 1 and covering the range from to 0.99. In §3.1, 
we create an equilibrium sequence for each core mass 
by slowly scanning over the binary separation, thereby 
identifying the separations of first contact, of the secular 
instability (if it exists), and of the Roche limit. In §3.2, 
separate dynamical calculations of various initial separa- 
tion r then allow us to test the dynamical stability of the 
contact configurations and to follow any mass transfer 
and the merger in unstable systems. 

3.1. Equilibrium Sequences 

Representative snapshots along the nic = 0.1 equilib- 
rium sequence are presented in Figure 3. The structure 
of these solutions is shown both in projection onto the 
orbital plane (the xy plane) and in terms of the effective 
potential $e- The thick solid curves in Figure 3 are the 
surfaces of constant effective potential $e that mark the 
inner and outer Lagrangian surfaces. For fixed x, $e is 
minimum on the binary axis {y = z = 0), and this min- 
imum value is given as a dashed curve in Figure 3. In 
hydrostatic equilibrium, the fluid flUs up to a constant 

level $e*' that is independent of x. 

Referring to Figure 3 we see that at the initial separa- 
tion in our scan, r — 3.30, the system is tidally bulged 
and the binary is detached: the fluid does not extend out 
to the inner Lagrangian surface (also known as the Roche 
lobe). At a separation r = 2.71, the binary stars fill the 
inner Lagrangian surface and make first contact through 
the LI point, located at the origin. Once the separation 
decreases to r = 2.42, the binary reaches the secular in- 
stability limit, as marked by a minimum in the energy 
E and angular momentum J along this equilibrium se- 
quence (see below). The critical separation r = 2.37 
at which mass transfer commences is identified with dy- 
namical calculations described in §3.2; in this and other 
scans through equilibrium configurations, however, we 
suppress the mass transfer by enforcing symmetry in par- 
ticle properties (see §2.3.1). At a separation r = 2.24, the 
fluid extends out to the outer Lagrangian surface, mark- 
ing the Roche limit. At even smaller separations, for 
example at the separation r = 2.22 shown in the figure, 
no equilibrium configurations exist and the stars shed 
mass through the outer Lagrangian surface near the L2 
point. The variation of critical effective potentials and 
the degree of contact rj along this rric = 0.1 equilibrium 
sequence is illustrated in Figure 4. For future reference, 
we note the approximately linear dependence of rj on the 
separation r. 

Prom Figure 5, we note that as rric increases toward 1, 
our results approach the Keplerian solution of two orbit- 
ing point masses. As expected, deviations from the point 
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Fig. 3. — Sequence of binary equilibrium configurations for two 
identical condensed polytropes of core mass nic = 0.1. Projections 
onto the orbital plane (the xy plane) are shown at six different 
binary separations for those SPH particles with \z\ < 0.06. The 
thick solid curves represent two surfaces of constant effective po- 
tential "I>e (see eq. [4]): namely, the inner and outer Lagrangian 
surfaces passing through the points Ll and L2. Shown beneath 
each configuration are corresponding projections onto the (a;, <J>e) 
plane for the same particles. The dashed curves give the vajriation 
of $e along the binary axis {y = z = 0). Contact configurations 
are obtained when the binary separation r < 2.71 (in units where 
an isolated binary component has radius R = 1). For r < 2.24, 
mass shedding through the outer Lagrangian points occurs, and 
no equilibrium configuration exists. 

mass result increase as a given binary becomes more 
deeply in contact or as we consider a binary associated 
with a smaller core mass. From the bottom frame of Fig- 
ure 5, we note that smaller core masses (corresponding 
to stars with less centrally concentrated density profiles) 
have a smaller orbital period at any given separation. For 
mc < 1, tidal interactions between the two stars make 
the effective potential stronger than 1 /r and shorten the 
rotation period compared to a point mass system. For 
rric = 0.1, for example, the deviation of the orbital pe- 
riod from the point mass result is approximately 1% at 
first contact, 2% at the secular instability limit, and 3% 
at the Roche limit. 

As in Rasio & Shapiro (1995), we determine the secu- 
lar stability limit along the equilibrium sequence by lo- 
cating where both the total energy E and total angular 
momentum J are minimum in curvets such as those of 
Figure 5. Our numerical results provide an accurate de- 
termination of this point for a close binary system, as the 
separate minima in E and J coincide to high numerical 
accuracy. This is in accord with the general property 
that dE = Qorh dJ along any sequence of uniformly ro- 
tating fluid equilibria (Ostriker & Gunn 1969). 

For the rric = binary, secular instability occurs soon 
after contact along this sequence and therefore stable, 
long-lived equilibrium configtirations can exist only in 
shallow contact, r] < 0.2. In contrast, the sequences with 
non-zero core masses permit much deeper contact before 
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Fig. 4. — Variation of critical effective potentials along the equi- 
librium sequence presented in Fig. 3. Values of the effective poten- 
tial at the outer Lagrangian surface (solid curve o), the inner La- 
grangian surface (solid curve i), and the fluid surface (long dashed 
curve s) are showrn as a function of binary separation r in the top 
frame. The degree of contact 7] (eq. [5]) is shown in the bottom 
frame as the long dashed curve. The short dashed curves give the 
positions of first contax;t {rj = 0) and of the Roche limit (77 = 1). 



the secular instability is reached. For example, a binary 
with core masses of 0.125 does not reach the secular in- 
stability until nearly r] = 0.9. For core masses rUc > 0.15, 
the stars will reach the Roche limit before the secular in- 
stability limit. Tables 2, 3, and 4 respectively present 
system properties at first contact [r] = 0), the secular 
instability limit (where E and J are minima), and at the 
Roche limit (77 = 1) for our sequences of various core 
mass rric. Additional runs at varying resolution indicate 
that the results in our tables have converged to within 
1% (e.g. see Fig. 6). 

These results for critical sciparations arc summarized 
in Figure 7. Due to tidal effects, the volume of each star 
is typically 1-2% larger at first contact than it is for that 
star in isolation. The separation Vfc at first contact is 
only weakly dependent on rric, being within 2% of 2.7 
for any core mass. For comparison, we note that the 
standard simple treatment of twin binaries would imply 
a first contact separation of 1/0.3799 = 2.63, where the 
factor 0.3799 comes from numerical integration of Roche 
lobe volumes around identical point masses (e.g., Eggle- 
ton 1983) and any change in the volumes of the stars due 
to tidal effects is neglected. We see therefore that finite 
size effects act to increase the separation of first contact. 

All three of the critical separations considered (first 
contact, secular instability, and Roche limit) tend to de- 
crease as the core mass increases. It is straightforward to 
find fitted formulas for the critical separations that are 
accurate to within ^ 1% for any core mass rric- for first 
contact 

rfc w 2.66 + 0.08(1 - mc)^ (6) 



Fig. 5. Variation of the system energy E (relative to the total 
self-energy Ex of the binary components at inflnity), angular mo- 
mentum ,/, and orbital period P along the equilibrium sequence 
for twin binaries with nic = 0, 0.1, 0.25, 0.5, and 1. In the E and 
,/ frames, higher curves correspond to smaller core mass, while in 
the P frame higher curves correspond to larger core mass. The 
orbital period P is normalized to the analytic point mass result 
^kepier = 2^^^TTr^^^ . The curves are from SPH scans of the equi- 
librium sequence, except for the rUc = 1 curve which is the analytic 
result for two point masses. The rric = and 0.1 curves exhibit a 
minimum in E and J, marking the position of the secular instabil- 
ity limit. The curves from the SPH scans terminate at the Roche 
limit, where mass shedding through the outer Lagrangian point 
commences. Critical points along our SPH scans are presented in 
Tables 2, 3, and 4. The individual data points (open circles) in 
this figure result from relaxing for an additional 200 time units at 
the given separation r. The agreement between these points and 
their corresponding scan helps to confirm that our scans axe indeed 
producing equilibrium sequences. 



TABLE 2 

FIRST CONTACT ALONG THE EQUILIBRIUM 
SEQUENCES OF TWIN BINARIES'' 





r 


P 


E- 


' Eoc 


J 


0.000 


2.75 


20.1 


-0 


154 


1.33 


0.050 


2.73 


19.8 


-0 


159 


1.31 


0.100 


2.72 


19.7 


-0 


162 


1.29 


0.125 


2.72 


19.8 


-0 


163 


1.29 


0.150 


2.71 


19.7 


-0 


165 


1.28 


0.175 


2.69 


19.5 


-0 


167 


1.27 


0.200 


2.69 


19.5 


-0 


168 


1.27 


0.250 


2.69 


19.5 


-0 


170 


1.25 


0.300 


2.70 


19.6 


-0 


171 


1.25 


0.400 


2.68 


19.4 


-0 


175 


1.23 


0.500 


2.66 


19.3 


-0 


178 


1.21 


0.600 


2.67 


19.3 


-0 


180 


1.20 


0.700 


2.66 


19.2 


-0 


183 


1.18 


0.800 


2.67 


19.4 


-0 


184 


1.17 


0.900 


2.67 


19.4 


-0 


186 


1.16 


0.990 


2.70 


19.7 


-0 


185 


1.16 



''Units arc defined such that G = M = i? = 1, rric is the core 
mass of each component, r is the binary separation, rj is the degree 
of contact (eq. [5]), P is the orbital period, and E — Eoo and J 
are the orbital energy and angular momentum, respectively; the 
energy E^ is the total equilibrium energy at infinite separation 
(that is, twice the energy of a single component in isolation). 
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TABLE 3 

SECULAR INSTABILITY ALONG THE EQUILIBRIUM 
SEQUENCES OF TWIN BINARIES'^ 



nic 


r 


V 


P 


E — Eoo 


J 


0.000 


2.67 


0.17 


19.1 


-0.155 


1.33 


0.050 


2.55 


0.40 


17.7 


-0.161 


1.30 


0.100 


2.42 


0.64 


16.4 


-0.167 


1.28 


0.125 


2.29 


0.86 


15.1 


-0.171 


1.26 


0.150 


2.22 


1.0 


14..3 


-0.175 


1.25 



''In twin binary sequences with core masses rric > 0.15, the Roche 
limit is reached before the secular limit. Units and column headings 
are as in Table 2, footnote a; the degree of contax;t r] is defined by 
eq. 5. 



TABLE 4 

ROCHE LIMIT^^ ALONG THE EQUILIBRIUM SEQUENCES 
OF TWIN BINARIES'^ 



rric 


r 


P 


E — A'oc 


J 


0.000 


2.34 


15.2 


-0.135 


1.38 


0.050 


2.28 


14.7 


-0.153 


1.32 


0.100 


2.24 


14.5 


-0.166 


1.28 


0.125 


2.22 


14.3 


-0.171 


1.27 


0.150 


2.21 


14.3 


-0.175 


1.25 


0.175 


2.20 


14.2 


-0.179 


1.24 


0.200 


2.19 


14.1 


-0.183 


1.22 


0.250 


2.17 


14.0 


-0.189 


1.20 


0.300 


2.16 


13.9 


-0.195 


1.18 


0.400 


2.14 


13.8 


-0.204 


1.15 


0.500 


2.13 


13.7 


-0.212 


1.12 


0.600 


2.12 


13.6 


-0.219 


1.10 


0.700 


2.12 


13.6 


-0.224 


1.08 


0.800 


2.10 


13.5 


-0.231 


1.05 


0.900 


2.10 


13.5 


-0.235 


1.04 


0.990 


2.12 


13.7 


-0.236 


1.03 



° Defined as the equilibrium configuration with the minimum bi- 
nary separation. 

"^Unit and column definitions are identical to those in Table 2, 
footnote a. 
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Fig. 6. — Critical separation r, energy E — Eao, angular mo- 
mentum J, and orbital period P at first contact (squares), secular 
instability (triangles), and the Roche limit (circles) versus total 
particle number for several rric = 0.1 equilibrium scans. Note 
the convergence of results for N > 10'* . Most of the binary calcu- 
lations in this paper employ AT ss 4 x 10'* particles. 
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Fig. 7. Separation at first contact (squares), the onset of secular 
instability (stars), and the Roche limit (circles) versus core mass 
rric- The data points arc determined from SPH scans of equilibrium 
sequences, as reported in Tables 2, 3, and 4. The lines are simply 
to help guide the eye. 



and for the secular instability limit 
rsen « 2.69 - Snir 



(7) 



We note that the 1 — rric in equation (6) equals the enve- 
lope mass (in units where the total stellar mass M — 1). 
We give our fit for the Roche limit separations in the 
next subsection, where we can determine these data with 
slightly better accuracy. We have not fit for the slight 
increase in the first contact data as the core mass is in- 
creased from rric = 0.9 to 0.99, as this feature is a numer- 
ical artifact due to our nic = 0.99 single star equilibrium 
model settling to a radius a few percent larger than 1. 

The critical core mass rric ~ 0.15, for which the secular 
instability and Roche limits coincide, can be determined 
graphically from Figure 7 by extrapolating the line con- 
necting the secular instability data down to the Roche 
limit curve. This intersection is important because it im- 
plies that all of our equal mass binaries with nic ^0.15 
can stably exist in deep contact, at separations all the 
way down to the Roche limit. Thus, essentially all twin 
red giant binaries will coalesce only due to mass shed- 
ding through the outer Lagrangian points. In contrast, 
twin binaries with mc ^ 0.15, corresponding primarily 
to main-sequence stars and subgiants, reach the secular 
instability limit at a larger separation than that of the 
Roche limit. As we will see in the next subsection, the 
secular instability limit is usually accompanied by dy- 
namical mass transfer at the same or a slightly smaller 
separation. Therefore, main sequence and subgiant twin 
binaries, as contrasted to most red giant twins, will start 
coalescing (a) when in more shallow contact and (b) 
through mass transfer across the inner Lagrangian point. 

3.2. Dynamical Integrations 
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Fig. 8. — Results of dynamical integrations: stable binaries (open 
squares), dynamically unstable binaries (asterisks), and configura- 
tions with no equilibrium (filled squares) . The curves represent the 
critical separations as determined by scanning runs, as in Fig. 7. 
The agreement between the results of the equilibrium sequences 
and dynamical calculations is excellent. 



Wc now study the stability of binary configurations 
with fully dynamical SPH integrations (see §2.4 for de- 
tails of the setup). Figure 8 summarizes the results of 
nearly 100 dynamical simulations with various rric and 
initial r values.^ We find the Roche limit to be very 
nearly at the separations determined from the equilib- 
rium scans, and the additional relaxation that we per- 
form before beginning a dynamical calculation allows us 
to determine these separations even more accurately. In 
addition, we find that most systems that are secularly 
unstable arc also dynamically unstable to mass transfer 
and then merger, as discussed below. 

The time evolution of the separation of nic = 0.05 
twin binaries is illustrated in Figure 9 for several dif- 
ferent initial separations. This figure also indicates the 
secular instability limit r^cc = 2.547, as determined by 
the energy and angular momentum minima in the equi- 
librium sequence of these binaries (see §3.1). Systems 
with separations r > Tgec are clearly dynamically stable, 
while those with r < rscc are dynamically unstable. That 
is, the secular and dynamical stability limits coincide or 
very nearly coincide, at least at this core mass. 

The bottom plot of Figure 9 shows the dynamical evo- 
lution of two cases that straddle the instability limit. For 
the r = 2.56 case, an epicyclic period of 350 time units 
is clearly evident, much larger than the orbital period 
of 17.9 time units. The large difference in these periods 
is an indication of how close the system is to an insta- 
bility limit (Rasio & Shapiro 1994). Indeed, if r were 
precisely at the dynamical stability limit, the period of 
small epicyclic oscillations would formally be infinite. 

^ Visualizations of selected simulations are available at 
http: / /webpub. allegheny.edu / employee/j /j alombar / movies / . 
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Fig. 9. — Top plot: Separation r versus time in several different 
dynamical SPH calculations for the core mass rric = 0.05 twin 
binaries. The initial separations are r = 2.64, 2.62, 2.56, 2.54, 
2.40, 2.30 and 2.28. Bottom plot: Zoomed in view of the r = 2.56 
and 2.54 cases, which straddle the instability limit. In both plots, 
the horizontal dotted line represents the secular instability limit at 
r = 2.547, as identified by an equilibrium scan. 



Figure 10 presents projected particle positions (top) 
and column density plots (bottom) at six different times 
in the nic = 0.05, r = 2.54 dynamical calculation, a 
case just inside the instability limit. Colors in the parti- 
cle plot are used to indicate from which component the 
particles originated. The coordinate system used here ro- 
tates counterclockwise with a period of 17.64 time units, 
which equals the orbital period of this binary at early 
times, before significant mass transfer has occurred. The 
instability initially manifests itself in the form of a nar- 
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row arm of gas that begins in the outer layers of one star, 
gradually flows across the neck surrounding the inner La- 
grangian point, and then creeps around the other star 
(see the t = 746 and 771 particle plot frames). The mass 
transfer drives the binary components closer, trigger- 
ing the excretion of mass through the outer Lagrangian 
points {t = 771) and accelerating the inspiral of the two 
cores. By t — 784, the mass transfered from one star 
has completely engulfed the other. At later times, the 
merger product approaches a rapidly rotating, axisym- 
metric configuration centered on the two cores orbiting 
in a tight binary (see the t = 792 and 961 frames). 

Figure 11 shows the evolution of binary separation for 
rric = 0.2 cases. Recall from §3.1 that there is no sec- 
ular instability limit at this core mass. Instead, stable 
equilibrium models exist all the way to the Roche limit. 
A binary at an initial separation r = 2.22 (or larger) 
orbits stably. In contrast, a binary at r = 2.20 grad- 
ually loses mass through the outer Lagrangian points, 
triggering a stage of rapid coalescence. The period of 
mass loss persists for several orbital periods: each of the 
oscillations superposed on the decreasing r = 2.20 curve 
in the bottom plot of Fig. 11 corresponds to one orbital 
period. These dynamical calculations indicate that the 
Roche limit for a rric = 0.2 twin binary indeed occurs 
near r = 2.2, in excellent agreement with the r = 2.19 
critical value estimated from the equilibrium scan. 

Figure 12 presents both particle positions and column 
density plots at six different times in the nic = 0.2, 
r = 2.20 dynamical calculation, a case just inside the 
Roche limit. The coordinate system used here rotates 
counterclockwise with a period of 14.22 time units, which 
equals the orbital period of this binary at early times. 
Gas is excreted almost immediately, with each parcel 
of gas carrying a specific angular momentum essentially 
equal to that of the outer Lagrangian points. In contrast 
to mergers with rric ^ 0.15, the excreted gas originates 
equally from both binary components and flows past the 
outer Lagrangian points symmetrically. As the outer La- 
grangian points are the outermost positions at which gas 
can be in static equilibrium, they are also the positions of 
largest possible speciflc angular momentum in rigidly ro- 
tating equilibrium twin binaries. Consequently, the mass 
loss necessarily decreases the average angular momentum 
per unit mass of the gas remaining within the outer La- 
grangian surface, causing the binary components, along 
with their cores, to inspiral: see the appendix of Webbink 
(1976) for a rigorous analysis of the angular momentum 
budget during mass excretion. As the components get 
closer, the excretion rate increases, and in addition the 
resulting arms become more tightly wound (compare the 
t = 202 frame to later ones). By t = 233, the central re- 
gions of the binary components have effectively merged. 
At later times, the merger product approaches a rapidly 
rotating, axisymmetrie configuration (see Fig. 13). 

Figure 14 shows energies versus time for the rric = 0.2, 
r = 2.20 calculation. The kinetic energy T gradually 
increases as the binary components inspiral, until the 
cores approach closely at t « 230. The ensuing shocks 
cause the gas to expand, causing an overall decrease in 
the internal energy U and increase in the gravitational 
potential energy W. The rapid variations in T and W 
at late times are due to the eccentric orbit of the central 
double core. The total energy E is well conserved in this 



simulation, varying by only 0.4% from its minimum to 
maximum values over the entire time interval shown in 
Figure 14. Energy conservation in other runs is typically 
at least this good and often even much better. In our 
simulations, most of the small non-conservation in energy 
occurs at late times once the core particles have entered 
a tight orbit. 

In all our merger simulations of twin binaries with frac- 
tional core masses of 0.15 or larger, each star loses mass 
through an outer Lagrangian point. Most of this mass 
ultimately ends up in a circumbinary envelope gravita- 
tionally bound to the central binary of cores. The mass 
ejected varies from ^ 0.5% (for m,c = 0.15) up to ^ 7% 
(for rric = 0.8) of the total system mass. A trend ev- 
ident from the simulations is that more massive cores 
ultimately remain more widely separated, after inspiral- 
ing, than less massive cores. The top frame of Figure 15 
provides a closer look at how the separation of the cores 
evolves in several dynamical simulations that begin just 
inside the Roche limit. For very large fractional core 
masses (rric ^ 0.9), the gas is simply not massive enough 
to affect significantly the dynamics of the cores: although 
they inspiral, they stay separated at distances on the or- 
der of the initial stellar radius. For moderately large core 
masses (0.5 < rric < 0.9), the cores inspiral to a fraction 
of a stellar radius, although the process halts at separa- 
tions large enough still to be resolved by our simulations. 
For core masses in the range 0.15 < mc < 0.5, which cor- 
responds to most red giants in nature, the cores rapidly 
inspiral to separations less than 0.1. Although our code 
has no mechanism for merging the cores, we do not ex- 
pect such an effect to be relevant here: the size of a core 
relative to the stellar radius is typically only ~ lO"** or 
less for a giant, so that a tremendous amount of angular 
momentum would have to be removed from the double 
core before they could merge. 

The bottom frame of Figure 15 concerns fractional core 
masses less than 0.15, namely those cases that reach the 
secular instability limit before the Roche limit. The 
cores in these merger simulations inspiral to a separa- 
tion < 0.01, considerably less than the spatial resolution. 
Whether or not the cores would merge in such circum- 
stances will likely depend on the details of the parent star 
structure, with simulations that resolve the cores neces- 
sary to address the issue fully. We find that less than 
0.5% of the system mass is ejected whenever the merger 
is triggered by mass transfer. 

Figure 16 helps to quantify the entropy evolution dur- 
ing mergers by plotting the mass-average < InA > over 
time for several cases. We note that, for our poly tropic 
equation of state and gas of uniform composition, the 
speciflc entropy of a parcel of gas is proportional to In A 
plus a constant. From the curves of Figure 16, we see 
that the change in entropy per unit mass tends to be 
larger, and develops on a longer timescale, for cases in- 
volving larger core masses. We also note that the entropy 
is still gradually increasing even at the end of our simula- 
tions, due to the influence of the central binary embedded 
within the system. 

Like the critical separations for first contact and sec- 
ular instability, the Roche limit separation tends to de- 
crease as the core mass increases. A fitted formula con- 
sistent with our dynamical integrations to within ~ 1% 
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Fig. 10. — Top frames: Projection of SPH particles onto the orbital plane at six different times in the merger of a rric = 0.05 binary with 
initial separation r = 2.54. Particles are colored according to the star in which they originated. Bottom frames: Column density plots at 
the same six moments, with the asterisks representing the positions of the compact cores. 



for any core mass rric is 

rRocho« 2.11 + 0.25(1 -mj^ 



(8) 



Because the degree of contact rj varies nearly linearly 
with separation r (for two examples, see Fig. 5 of this pa- 
per and Fig. 2 of Rasio & Shapiro (1995)), this formula, 
along with others from §3.1, allows us also to estimate 
the degree of contact rj at the secular instability limit: 

'7soc ~ {rfc - ?'scc)/(?'fc - r-Rochc)- 

4. DISCUSSION AND FUTURE WORK 

4.1. Summary of Main Results 

We have determined equilibrium sequences and per- 
formed dynamical calculations of twin binaries, focusing 
primarily on configurations in which the stars are in con- 
tact. Our equilibrium sequences of §3.1 allow us to de- 
termine the binary separation at first contact and at the 
innermost stable circular orbit as a function of the frac- 



tional core mass rric. For rric ^ 0.15, the innermost stable 
orbit occurs at the secular instability limit (marked by 
a minimum in energy and angular momentum along the 
equilibrium sequence); for rUc > 0.15, the innermost sta- 
ble orbit occurs at the Roche limit (defined as the min- 
imum separation for which an equilibrium configuration 
exists). Our dynamical calculations of §3.2 confirm these 
critical separations and also reveal how the components 
inspiral once a binary passes the innermost stable orbit. 

Figure 17 summarizes graphically our most basic re- 
sults. Recall that the separation r on the vertical axis 
is scaled to the unperturbed stellar radius R, while the 
core mass rric on the horizontal mass is scaled to the stel- 
lar mass M. Thus, as the components in a twin binary 
expand and gradually increase their core masses due to 
stellar evolution, the corresponding position in the pa- 
rameter space of Figure 17 will shift down and slightly 
to the right. When this position drops below the top 
curve, which marks first contact, the binary enters the 



Twin Binaries 



13 



2.5 - 



2 - 



1.5 - 




1 - 



100 200 300 

t 




2.1 - 



100 200 300 400 500 

t 

Fig. 11. — Like Fig. 9, but for twin binaries with core mass 
rric = 0.2. At this core mass, there is no secular instability limit: 
instead, the system become unstable only once it reaches the Roche 
limit at r ^ 2.2. The initial separations shown in the top plot are 
r = 2.60, 2.40, 2.22, 2.20, while a zoomed in view of the r = 2.22 
and 2.20 curves are shown in the bottom plot. 

stable contact phase. When the position drops into ei- 
ther the unstable or no equilibrium portions of parameter 
space, the components merge. 

The curves shown in Figure 17 are given by equations 
(6), (7), and (8); population modelers can use these fitted 
formulae in treatments of twin binaries. Consider, for 
example, such a binary with a given orbital separation 
a. The stellar evolution of each component gives the 
time dependence of the stellar radius R, the stellar mass 
M, and the core mass Mg. The dimensionless separation 
r = a/R and the fractional core mass rric = Mc/M are 



thus known functions of time, and the evolutionary track 
can be placed in the theoretical r versus nic diagram 
(Figure 17) to determine the final fate of the system. 

We note that the volume of parameter space where real 
binaries would ultimately end up in the stable contact 
region, without crossing the instability limit or the Roche 
limit, is small. Such a situation would require a fine 
tuning of the initial semimajor axis a. For example, a 
star with an initial mass of 8Mq will expand to i? « 
37Oi?0 and reach a fractional core mass of nic ~ 0.2 
as it ascends the asymptotic giant branch, according to 
calculations by the TWIN stellar evolution code. For 
this core mass, rfc ~ 2.7 and rRochc ~ 2.2. Thus, a twin 
binary composed of such stars will remain detached if a > 
rfcR « IOOORq and will ultimately surpass the Roche 
limit if a < r^ocheR » SOORq. Only if SOOi?© < a < 
lOOOiiQ will the binary reach the contact phase without 
the cores also inspiraling to form a tight binary. 

The dynamic simulations of §3.2 always start from 
a symmetric equilibrium configuration, with the binary 
components being hydrostatic in the corotating frame. 
In our simulations with core masses of nic = 0.125 and 
less, we find a dynamical instability to exist at or slightly 
inside the secular instability limit. This dynamical in- 
stability is a global instability of the equilibrium state, 
triggered by small numerical noise and characterized by 
a growing asymmetric mode. Binaries that come in- 
side this instability limit first transfer mass gradually 
from one component to the other and eventually coa- 
lesce quickly as mass is lost through the outer Lagrangian 
points. The cores are left in a tight binary surrounded 
by a circumbinary envelope. 

4.2. Comparison with Other Works 

The merger of our small core mass twin binaries pro- 
ceeds in a fashion qualitatively similar to that of the 
Q1.3 model of D'Souza et al. (2006). In that model, 
the binary is composed of purely polytropic components 
{nic = 0) with mass ratio (donor to accretor) q = 1.3. In 
both our simulations and theirs, the dynamical instabil- 
ity manifests itself as a gradually developing mass trans- 
fer flow, followed by excretion of gas through the outer 
Lagrangian points and merger of the stellar envelopes. 
One important difference, however, is that the instabil- 
ity in our rric = case does not develop until the stars 
have reached a degree of contact ry w 0.17, whereas the 
instability is present in the Q1.3 model while the binary 
is still semidetached. This difference highlights the stabi- 
lizing influence of the common envelope in twin binaries, 
even though the instability still exists even for g = 1. 

We can directly compare our results for the limiting 
case rUc = with those of Rasio & Shapiro (1995). 
The agreement in the first contact, secular instability, 
and Roche limit separations is excellent (better than 
1%). The computational resources of the time, how- 
ever, limited Rasio & Shapiro (1995) to follow up to 
only ~ 3 orbits, so that they were unable to identify 
weak mass transfer events, that is, events that develop 
gradually over many dynamical timescales (Paczyinski 
& Sienkiewicz 1972). As a result, they determined the 
dynamical instability limit to be at r w 2.45, well inside 
the secular instability limit. In contrast, by following the 
dynamical evolution of nic = twin binaries for up to 
~ 100 orbits, we find that the dynamical instability limit 
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Fig. 12. — Like Fig. 10, but for the rUc = 0.2 dynamical calculation starting at r = 2.20, just inside the Roche limit. 



actually coincides, or nearly coincides, with the secular 
instability limit at r « 2.67. 

4.3. Relevance to Binary Neutron Stars 

We find that twin binaries with rUc > 0.15 exist stably 
at separations all the way down to the Roche limit, where 
mass is then excreted symmetrically through the outer 
Lagrangian points. This excretion carries away angular 
momentum and causes the stars, along with their cores, 
to inspiral on a dynamical timescale. For core masses 
nT-c ^ 0.9, the cores inspiral to a final separation that 
is a fraction of the stellar radius. Indeed for mc < 0.5, 
which corresponds to most giant stars in nature, the final 
separation of the cores is less than one-tenth the stellar 
radius. Thus, we are left with two cores in a tight binary 
surrounded by the combined gaseous envelopes of the 
original binary, the precursor for double neutron stars 
proposed in the Brown (1995) scenario. 

As the gravitational forces of the cores are softened at 
distances less than ~ 0.1, our dynamic simulations that 



lead to cores in a tight binary can provide only an upper 
limit on the separation at the end of their inspiral. Our 
simulations of the rUc = 0.15 case, for example, place 
this upper limit at ^ 0.03 times the stellar radius (see 
Fig. 15). Thus, a binary composed of twin M — IOMq, 
R — 2OOi?0, Mc — 1.5M0 red giants would have their 
cores spiral to a separation of less than 0.03i? ~ 6Rq. 
The gradual transfer of energy to the circumbinary en- 
velope could easily decrease the separation of the cores 
by a factor of '--^ 2 further. As gravitational radiation 
alone can bring two 1.5M0 point masses separated by 
up to ~ 5i?0 to contact in less than a Hubble time, such 
systems could have evolved to become arbitrarily tight 
by the present time. We therefore believe that such dou- 
ble cores are indeed excellent candidates for binary NS 
progenitors, as proposed in the Brown scenario. 

4.4. Relevance to Planetary Nebulae 

Given the relatively short timescale covered by hydro- 
dynamic simulations such as ours, the circumbinary en- 
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Fig. 14. — Energies versus time t for the simulation presented in 
Fig. 12. From the bottom curve to the top one, we show gravi- 
tational energy W, total energy E, kinetic energy T, and internal 
energy U. 
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of time t for binaries that begin just inside the Roche limit (top 
frame) and those that begin just inside the secular limit (bottom 
frame). Each curve is labeled by the factional core mass rric. The 
horizontal dashed line indicates the minimum separation for which 
the gravitational interaction of the two cores is unsoftened in our 
simulations. 



velopes at the end of our dynamical simulations are still 
optically thick. Nevertheless, it is natural to think of 
the final states of our merger calculations as a type of 
proto-PNe, specifically, as the immediate precursors of 
PNe with equal mass central binaries. Future work could 
use the end state of hydrodynamic calculations as initial 
conditions in models of PN evolution, with particular at- 
tention paid to any transition to the optically thin regime 
(revealing the central binary) and to the morphology of 
the gas distribution. 



To examine whether the dynamical calculations of this 
paper yield final states that are consistent with the obser- 
vations of the g « 1 double degenerate binaries, we con- 
sider a twin binary with 3Mq components. We use the 
TWIN stellar evolution code to evolve each star in isola- 
tion. The calculation shows that the star expands to only 
R « 30i?o as it ascends the giant branch and to nearly 
5OOi?0 on the asymptotic giant branch. Therefore, for 
any reasonable distribution of initial orbital separations, 
twin binaries made of such stars would much more often 
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Fig. 16. — A measure of entropy evolution due to shock heating 
during several mergers. The natural logarithm of tlie entropic vari- 
able A is averaged by mass over the gas of the system and plotted 
versus a shifted time coordinate At, chosen such that the inflection 
point in each curve (when the entropy is increasing most rapidly) 
occurs at At = 0. Each curve is labeled by the factional core mass 
mc- as rUc increases, so do both the initial and final values of 
< In A >. For each core mass, we present our simulation of largest 
initial separation that results in a merger, corresponding in Fig. 8 
to the uppermost asterisk (for a given mc < 0.15) or filled square 
(for a given rric > 0.15). Note that the initial value of the entropic 
variable A equals the K in equation (2), as calculated from the 
appropriate Eq in Table 1. The entropic variable A is in units of 



Fig. 17. — The parameter space of twin binaries. Here and 
throughout the paper, the core mass rric and binary separation 
r arc given in units of the total mass and radius, respectively, of 
an isolated binary component. The top curve, separating detached 
and contact binaries, marks configurations of first contact. The 
middle curve, spanning < rtic < 0.15, is the secular instability 
limit that separates stable and unstable conta<:t binaries: in this 
work, we find that most secularly unstable systems are also dy- 
namically unstable to mass transfer across the inner Lagrangian 
point. The bottom curve, separating contact binaries from config- 
urations with no equilibrium, represents the Roche limit. Binaries 
that are unstable or that cannot exist in equilibrium have their 
components inspiral and merge, a process we follow with dynami- 
cal calculations. 



come into contact while on the asymptotic giant branch 
than when on the red giant branch, consistent with the 
fact that most or all of the well-observed double degen- 
erate stars in PNe seem to have masses too large to be 
He white dwarfs. 

For the sake of discussion, consider such a twin binary 
that reaches the Roche limit when R — 2OOi?0. Ac- 
cording to the stellar evolution calculation, an isolated 
star at that time would have a carbon-oxygen core of 
mass Mc = O.STM© and, accounting for stellar winds, 
a total mass M = 2.8Mq. At this fractional core mass 
rric = 0.57/2.8 = 0.20, the dimcnsionless Roche separa- 
tion (from the results of this paper) is rRoche ~ 2.2, and 
thus the initial semimajor axis of this binary would have 
been a = TRoche-R ~ 44Oi?0. From the rric = 0.2 curve in 
the top frame of Figure 15, we see we can place an upper 
limit of ^ 0.05i? = lOi?0 on the final separation of the 
two cores, which, from Kepler's third law, corresponds to 
an upper limit on the orbital period of 3 days. If instead 
R — lOOi?0 at the time of merger, then a similar calcula- 
tion gives a core mass Mc = O.SSMq and an upper limit 
on the orbital period of 0.9 days. 

For comparison, the central binaries in NGC 6026, 
Abcll 41, and Hen 2-428 have comparable masses to these 
core masses and orbital periods in roughly the 0.2 to 0.5 
day range. We conclude that the results of our dynami- 
cal simulations, when applied to intermediate mass twin 
binaries, are consistent with observed characteristics of 



double degenerate central binaries in PNe. In future sim- 
ulations, a reduced gravitational softening between the 
two cores could allow for an even more precise determi- 
nation of their final orbital separation. In this way, ob- 
served orbital parameters of degenerate binaries in PNe 
could be more readily connected to the properties of the 
parent stars from which they may have originated. 

4.5. Additional Comparisons and Future Work 

For core masses rUc > 0.15, we find that a twin bi- 
nary can exist stably in deep contact, at separations 
all the way down to the Roche limit. In contrast, the 
semi-analytic condensed polytrope models of Hjellming 
& Webbink (1987) predict instead that twin binaries will 
experience sustained mass transfer once the components 
come in contact, provided only that rUc < 0.458 (a range 
that includes thc^ vast majority of giant stars) . The pri- 
mary oversimplification in the semi-analytic treatment 
appears to be the approximation that mass outside of 
the Roche lobe cannot help to contain the star within it. 
Our numerical calculations, however, model the common 
envelope that exists outside of the Roche lobe and that 
acts to suppress mass transfer. In addition, our fully 
three-dimensional calculations remove the point mass 
and spherical structure approximations implicit to the 
semi-analytic method. 

The models of Hjellming &; Webbink (1987) seem best 
suited to semidetached binaries, where there is no com- 
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mon envelope to complicate the dynamics of the mass 
flow. A comparison of such cases with our results is not 
possible, as our work is limited to binaries with identical 
components. Natural future work would include relax- 
ing this constraint so that binaries with mass ratio q^l 
can be studied and compared with semi-analytic models. 
Of particular interest to the binary neutron star problem 
would be cases in which the mass ratio deviated from one 
by only a few percent or less. 

The modeling of giants as F = 5/3 condensed poly- 
tropes is a common simplifying approximation, one that 
here allows our results to be scaled to binaries of any 
mass and length scales. We note, however, that radiation 
pressure can be the dominant contributor to the equation 
of state at some ages and at some locations within the 
envelopes of massive giants (M > IAMq, according to 
calculations with the TWIN stellar evolution code). For 
such cases, our treatment of the envelope as a constant 
entropy, T = 5/3 gas can be legitimately questioned. 
While the effects of employing more realistic stellar mod- 



els would be worthwhile to study in future work, we do 
not expect our results to change qualitatively. Regard- 
less of the equation of state, gas that flows past the outer 
Lagrangian points will still necessarily carry away a spe- 
cific angular momentum larger than the system average, 
forcing the remaining gas to configurations of smaller an- 
gular momentum per unit mass. We conclude that the 
inspiral of cores should be a common outcome whenever 
a real twin binary exceeds the Roche limit. 
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